We define both fuctions:

f[x1_,x2_,x3_,x4_]:={-2+x2 x4+2x3 x4+2(x4^2),-3(x1^2)-x1 x2+3x3+(x3^2)-2x4}

g[u_,v_]:={-u-3(u^2)+v-3u v,-1-(u^2)-2v+3u v+(v^2),3-3u-(u^2)+3v+3u v+(v^2),-2+2(u^2)-v-u v-2(v^2)}

If we apply the chain rule formula we have: (g∘f)’(p)=g’(f(p)f’(p). Therefore we need to compute: f’(p), f(p) and g’(f(p)).

As we are going to use matrices, we activate the option MatrixForm:

$Post=MatrixForm

We compute f’(p) by

D[f[x1,x2,x3,x4],{{x1,x2,x3,x4},1}]/.{x1->0,x2->1,x3->3,x4->3}

We compute now f(p)

f[0,1,3,3]

Finally we calculate g’(f(p))

D[g[u,v],{{u,v},1}]/.{u->37,v->12}

The Jacobian matrix of the composition is

(({{-259, -110}, {-38, 133}, {-41, 138}, {136, -86}})).(({{0, 3, 6, 19}, {-1, 0, 9, -2}}))

We end by computing its determinant

Det[(({{110, -777, -2544, -4701}, {-133, -114, 969, -988}, {-138, -123, 996, -1055}, {86, 408, 42, 2756}}))]

Therefore, the correct answer is option 4).

First of all we know a solution of the system given by point p=(x1,x2,x3,x4,u)=(5,4,-4,-5,-1). Therefore, for u=-1, we know that x1=5, x2=4, x3=-4, x4=-5 or, in other words if we solve x1, x2, x3, x4 in terms of variable u we will obtain the formulas x1(u), x2(u), x3(u), x4(u) for which we know the value for u=-1:
x1(-1)=5, x2(-1)=4, x3(-1)=-4, x4(-1)=-5

We first define the function F which groups all the equations of the system:

F[x1_,x2_,x3_,x4_,u_]:={-(u^2)x3+2u (x4^2)+42,3u x2 x3+(x2^2)x3+16,u x1 x2-u x1 x3+40,-3u x3+2(x2^2)x3+140}

We compute the Jacobian matrix at p:

D[F[x1,x2,x3,x4,u],{{x1,x2,x3,x4,u},1}]/.{x1->5,x2->4,x3->-4,x4->-5,u->-1}

The matrix corresponding to the variables we have to solve for Y=(x1,x2,x3,x4) is ((∂F)/(∂Y))(p):

(({{0, 0, -1, 20}, {0, -20, 4, 0}, {-8, -5, 5, 0}, {0, -64, 35, 0}}))

The matrix corresponding to the variables X=(u) is ((∂F)/(∂X))(p):

(({{42}, {-48}, {40}, {12}}))

The function with groups the formulas to compute x1, x2, x3, x4 in terms of u is G(u)=(x1(u),x2(u),x3(u),x4(u))  and the only value that we know for that functions is for p=(x1,x2,x3,x4,u)=(5,4,-4,-5,-1) so, for u=-1, we know that the values of x1,x2,x3,x4 are (x1,x2,x3,x4)=(5,4,-4,-5); that is to say G(-1)=(5,4,-4,-5)
The Jacobian matrix of G(u) at u=-1 is given by the implicit function theorem as
G(-1)= - (((((∂F)/(∂Y))(p)))^(-1))·((∂F)/(∂X))(p):

-Inverse[(({{0, 0, -1, 20}, {0, -20, 4, 0}, {-8, -5, 5, 0}, {0, -64, 35, 0}}))].(({{42}, {-48}, {40}, {12}}))

The tangent affine function at u=-1 of G(u) is
G(-1)+G’(-1)(u-(-1)):

(({{5}, {4}, {-4}, {-5}}))+(({{(205/74)}, {-(144/37)}, {-(276/37)}, {-(183/74)}}))(u-(-1));

In particular, in the exercise the ask to compute the approximation for x1,x2,x3,x4 when u=-0.8 so in the preceding expresion we substitute u with -0.8

(({{5}, {4}, {-4}, {-5}}))+(({{(205/74)}, {-(144/37)}, {-(276/37)}, {-(183/74)}}))(-0.8-(-1))

So the approximation is (({{x1}, {x2}, {x3}, {x4}}))=(({{5.554054054054054`}, {3.2216216216216216`}, {-5.491891891891892`}, {-5.494594594594594`}}))

So the correct answer is 4).
As the following exercises doesn’t need matrix calculus we disable the optin MatrixForm.

$Post=.

We have to compute the possible maxima/minima in two steps:
1- In the interior of the domain Int(D)≡((x^2)/9)+((y^2)/9)+((z^2)/25)<1
2- In the obundary of the domain ∂D≡((x^2)/9)+((y^2)/9)+((z^2)/25)=1

Step 1: interior of the domain

We can use the Jacobian matrix technique for the function

f[x_,y_,z_]:=-7+4x-(x^2)+4y-(y^2)-(z^2)

We solve the equation f’(x,y,z)=(0,0,0)

Solve[D[f[x,y,z],{{x,y,z},1}]=={0,0,0}]

We check that this point satisfy the equation of the domain

((x^2)/9)+((y^2)/9)+((z^2)/25)<1/.{x->2,y->2,z->0}

And we see in this way that this point is in the domain so we will consider it later.

Step 2: boundary of the domain

We have to compute the maxima/minima of f(x,y,z) with the constraint
((x^2)/9)+((y^2)/9)+((z^2)/25)=1. As we have an optimization problem with constraints we use the Lagrange function:

La[x_,y_,z_,λ_]:=f[x,y,z]-λ(((x^2)/9)+((y^2)/9)+((z^2)/25)-1)

For the Lagrange function we apply the Jacobian matrix technique: La’(x,y,z,λ)=(0,0,0,0):

N[Solve[D[La[x,y,z,λ],{{x,y,z,λ},1}]=={0,0,0,0}]]

We take only the solutions without complex numbers (in this case, all of them):

solutions={{x->-1.125`,y->-1.125`,z->-4.238956239453293`,λ->-25.`},{x->-1.125`,y->-1.125`,z->4.238956239453293`,λ->-25.`},{x->-2.1213203435596424`,y->-2.1213203435596424`,z->0.`,λ->-17.485281374238568`},{x->2.1213203435596424`,y->2.1213203435596424`,z->0.`,λ->-0.5147186257614291`}};

Let us substitute these points in the function to check where we obtain max/min value:

f[x,y,z]/.solutions

Also we have to consider the point computed in Step 1:

f[x,y,z]/.{x->2,y->2,z->0}

Therefore the maximum value of f(x,y,z) is 1 at point (2,2,0) and the correct options is 5).

We consider the point (x1,y1) in domain (Subscript[D, 1]) and (x2,y2) in domain (Subscript[D, 2]) and we have to compute the points (x1,y1) and (x2,y2) for which the square distance funcion is minimum:

d[x1_,y1_,x2_,y2_]:=(((x1-x2))^2)+(((y1-y2))^2)

For sure the minimum distance will be attained at the boundary (we can check it with a graphic representation of both domains) and therefore (x1,y1) has to satisfy the equation of domain (Subscript[D, 1]) and (x2,y2) the one of (Subscript[D, 2]) so we have two constraints and we need to use the corresponding Lagrange function

La[x1_,y1_,x2_,y2_,λ1_,λ2_]:=d[x1,y1,x2,y2]-λ1(1881+240x1+8(x1^2)+258y1+16x1 y1+9(y1^2)-1)-λ2(376-94x2+9(x2^2)+94y2-8x2 y2+7(y2^2)-100)

And we solve La’(x1,y1,x2,y2,λ1,λ2)=(0,0,0,0,0,0). We apply the instruction N to obtain the numeric expression instead of the exact one and also the instruction Chop to remove the small complex terms due to precission mistakes of Mathematica.

Chop[N[Solve[D[La[x1,y1,x2,y2,λ1,λ2],{{x1,y1,x2,y2,λ1,λ2},1}]=={0,0,0,0,0,0}]]]

We will consider only the solutions without complex terms:

solutions={{x1->-7.060165427871764`,x2->-0.8458794489611478`,y1->-8.067810670121908`,y2->-7.511701465203684`,λ1->6.069772086677446`,λ2->-0.252961750410661`},{x1->-4.956842876977498`,x2->-0.6871283615566625`,y1->-9.86694438409176`,y2->-8.2230367791704`,λ1->-3.028806643981546`,λ2->-0.21041360157777791`},{x1->-5.099093743542912`,x2->6.454367956098165`,y1->-9.624879674301056`,y2->-1.3408147693876118`,λ1->-5.232042149259568`,λ2->0.7022283582643081`},{x1->-7.025751872095876`,x2->6.620704149560221`,y1->-8.173034426637424`,y2->-1.622860568514477`,λ1->8.581109531368506`,λ2->0.715306303639486`}};

We compute de distance squares for all these points:

d[x1,y1,x2,y2]/.solutions

The largest smallest distance square is the second one which correspond to a distance

(Sqrt[20.93289425671366`])

And the corresponding point is also the second one

{x1->-4.956842876977498`,x2->-0.6871283615566625`,y1->-(Subscript[9.86694438409176`, \[Placeholder]]),y2->-8.2230367791704`,λ1->-3.028806643981546`,λ2->-0.21041360157777791`}

Therefore the smallest distance is obtained for (x1,y1)=(-4.956842876977498`,-(Subscript[9.86694438409176`, □])) in (Subscript[D, 1]) and (x2,y2)=(-0.6871283615566625`,-8.2230367791704`) in (Subscript[D, 2]). 
Therefore the correct option is 3) since the first coordinate of (x1,y1) which is -4.956842876977498 has digit 9 just after the dot symbol and it fits with 3).

In this case we have to compute the point (x,y) in domain (Subscript[D, 1])for which the distance to q=(-4,-8) is minimum. The square distance function is

d[x_,y_]:=(((x-(-4)))^2)+(((y-(-8)))^2)

We can check with a graphic representation that the minimun distance is reached at a point of the boundary and hece (x,y) hasto satisfy the equation of domain (Subscript[D, 1]):

6(x^2)+6x y+8(y^2)==1

So we have to compute a minimum with constraints and we need to use the Lagrange function:

La[x_,y_,λ_]:=d[x,y]-λ(6(x^2)+6x y+8(y^2)-1)

And for this function we use the Jaconbian matrix technique:

N[Solve[D[La[x,y,λ],{{x,y,λ},1}]=={0,0,0}]]

We consider only the real solutions:

solutions={{x->-0.07755978458509913`,y->-0.31924580349584664`,λ->-2.7562721630725737`},{x->0.06503299301341509`,y->0.3255022906516029`,λ->2.974331328666207`}};

And we compute for them the square distance:

d[x,y]/.solutions

The minimum square distance is 74.3795 which correspond to a distance

(Sqrt[74.37952227062026`])

And the point where we obtain the minimum distance to q is the first one

{x->-0.07755978458509913`,y->-0.31924580349584664`,λ->-2.7562721630725737`}

The distance 8.62436 fits with option 5).

To compute the diameter we have to computed the largest distance. We can check with a graphic that the maximum distance is obtained at the boundary of (Subscript[D, 1]). Therefore we have to compute the maximum square distance:

d[x_,y_,z_]:=(((x-1))^2)+(((y-0))^2)+(((z-(-1)))^2)

For a point (x,y,z) satisfying the equation of the boundary of (Subscript[D, 1]):

4(x^2)-4x y+8(y^2)+2x z-10y z+5(z^2)==1

Therefore we consider the Lagrange function

La[x_,y_,z_,λ_]:=d[x,y,z]-λ(4(x^2)-4x y+8(y^2)+2x z-10y z+5(z^2)-1)

And we use again the Jacobian matrix:

N[Solve[D[La[x,y,z,λ],{{x,y,z,λ},1}]=={0,0,0,0}]]

We take the real solutions:

solutions={{x->0.3773593397033856`,y->-0.1529820797206817`,z->-0.46221055008365425`,λ->-0.46012770920136725`},{x->-0.15359035087011186`,y->0.4190835469498634`,z->0.6755140800684764`,λ->1.484644718513582`}};

For these points the square distance is

d[x,y,z]/.solutions

The maximum square distance is the second one 4.31375 which correspond to a distance

(Sqrt[4.313749149452419`])

The corresponding point is the second one

{x->-0.15359035087011186`,y->0.4190835469498634`,z->0.6755140800684764`,λ->1.484644718513582`}

Which fits with option 3).

We have the points (x1,y1,z1) in (Subscript[D, 1]) and (x2,y2,z2) in (Subscript[D, 2]) and we have to compute the minimum square distance

d[x1_,y1_,z1_,x2_,y2_,z2_]:=(((x1-x2))^2)+(((y1-y2))^2)+(((z1-z2))^2)

By means of a graphics we can realize that the minimum distance is attained on the boundary so (x1,y1,z1) hasto satisfy the equation of domain (Subscript[D, 1]) and (x2,y2,z2) the one of (Subscript[D, 2]). So we have a problem with two constraints and the Lagrange function will have to Lagrange multipliers:

La[x1_,y1_,z1_,x2_,y2_,z2_,λ1_,λ2_]:=d[x1,y1,z1,x2,y2,z2]-λ1(9(x1^2)+7(y1^2)-6x1 z1+2y1 z1+2 (z1^2)-1)-λ2(8x2+3y2+14)

We solve the equation:

N[Solve[D[La[x1,y1,z1,x2,y2,z2,λ1,λ2],{{x1,y1,z1,x2,y2,z2,λ1,λ2},1}]=={0,0,0,0,0,0,0,0}]]

We discard any solution with complex numbers but also any solution with non numeric expressiones as the first ones (for instance, in the first one, the value for x1 is x1->0.001890359168241966` (-784.`+51.` z1-3.` (Sqrt[-11819.`-3948.` z1-769.` (z1^2)])) which is not a numeric expression as z1 appears in the formula). The only solutions that we will keep are ones with numeric non complexs values:

solutions={{x1->-0.4661906462318917`,x2->-1.586882087382779`,y1->-0.014721809881007127`,y2->-0.4349811003125899`,z1->-0.691925064407334`,z2->-0.691925064407334`,λ1->-0.5286428445511427`,λ2->-0.2801728602877218`},{x1->0.4661906462318917`,x2->-1.4816110633021524`,y1->0.014721809881007127`,y2->-0.7157038311942592`,z1->0.691925064407334`,z2->0.691925064407334`,λ1->0.9188001251193822`,λ2->-0.4869504273835111`}};

We compute the square distance for them:

d[x1,y1,z1,x2,y2,z2]/.solutions

The smallest one is the first one which correspond to a distance

(Sqrt[1.43256717746291`])

And the corresponding point is also the first one

{x1->-0.4661906462318917`,x2->-1.586882087382779`,y1->-0.014721809881007127`,y2->-0.4349811003125899`,z1->-0.691925064407334`,z2->-0.691925064407334`,λ1->-0.5286428445511427`,λ2->-0.2801728602877218`}

The point of (Subscript[D, 1]) with minimun distance will be (x1,y1,z1)=(-0.4661906462318917`,-0.014721809881007127`,-0.691925064407334`) which fits with option 4).